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Abstract 

The concept of nonlinear modes is applied in order to analyze the behavior of a 
model of woodwind reed instruments. Using a modal expansion of the impedance 
of the instrument, and by projecting the equation for the acoustic pressure on the 
normal modes of the air column, a system of second order ordinary differential 
equations is obtained. The equations are coupled through the nonlinear relation de- 
scribing the volume flow of air through the reed channel in response to the pressure 
difference across the reed. The system is treated using an amplitude-phase formu- 
lation for nonlinear modes, where the frequency and damping functions, as well as 
the invariant manifolds in the phase space, are unknowns to be determined. The 
formulation gives, without explicit integration of the underlying ordinary differen- 
tial equation, access to the transient, the limit cycle, its period and stability. The 
process is illustrated for a model reduced to three normal modes of the air column. 

Key words: nonlinear modes, model reduction, amplitude phase formulation, 
autonomous system, periodic oscillations, clarinet-like instruments 
PACS: 43.75.Pq, 43.25.Ts 



1 INTRODUCTION 

Musical wind instruments are interesting examples of nonlinear vibrating sys- 
tems. The process governing the formation of self-sustained oscillations is sur- 
prisingly complex. In short, a wind instrument consists of a resonator and a 
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generator [T]. The resonator is the air column inside the instrument, and is 
usually characterized by the linear wave equation. The generator, in turn, con- 
sists of some kind of pressure or flow controlled valve [2] , where the relationship 
between air flow and pressure is starkly nonlinear. Even a simplistic model of 
sound generation must include these nonlinear effects, simply because they are 
a prerequisite for the forming of a self-sustained oscillation from a continuous 
supply of air. One important reason why different wind instruments sound so 
different can be attributed to the nonlinearities, and how they interact with 
the vibrational modes of the air column at hand. The aim of this paper is 
to study how the sound production in a clarinet-like instrument can be ana- 
lyzed and simulated in the framework of nonlinear modes. A first reason for 
pursuing this subject is to evaluate the possibility to use nonlinear modes to 
derive models of reduced complexity, which are of interest for sound synthesis. 
Another goal is to identify important control parameters in a reduced model 
(functions of such entities as blowing pressure, pinching force and position of 
the player's lips on the reed etc) which can be regulated by a musician in an 
intuitive way, without a long period of training. 

Modal analysis is the natural tool for characterizing linear mechanical systems 
(in particular for numerical modelling, prediction and experimental character- 
ization). The extension of the modal theory to nonlinear mechanical systems 
appears for the first time under the name Nonlinear Normal Modes in the 
work of Rosenberg [3] for a system of n-masses interconnected by nonlinear 
springs. A NNM is defined as a family of periodic solutions of the equation 
of motions corresponding to simple curves in the configuration space. In this 
paper the name Nonlinear Modes (NM) will be used in place of NNM. Shaw 
and Pierre [I] extended the concept of NM in the context of phase space. This 
approach is geometric by nature and makes use of the theory of invariant man- 
ifolds for dynamical systems. A NM of an autonomous system is defined as a 
two-dimensional invariant manifold in the phase space, passing "through a sta- 
ble equilibrium point of the system and, at that point, it is tangent to a plan, 
which is an eigenspace of the system linearized about that equilibrium" [5]. In 
the invariant manifold, the modal dynamics reduces to a one-degree of freedom 
nonlinear oscillator. This definition is valid for dissipative mechanical systems. 
A nonlinear superposition technique is also proposed [5] and its validity is dis- 
cussed by Pellicano and Mastroddi [5]. The invariant manifold approach is 
also related to the methods based on the theory of the normal forms [8f9][T0] . 
where the invariant manifold and the modal dynamics equation of motion are 
extracted from the minimal representation. The construction of the NM for 
piecewise linear systems has been considered By Jiang, Pierre and Shaw [7]. 
A review paper by Vakakis [TT] and the book by Vakakis et al. [12] contain an 
almost complete report of the history of the subject. 

In this paper, a method devised by Bellizzi and Bouc [13] for computing 
two-dimensional invariant manifolds of dynamical systems is considered. This 
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method (recalled in section [3]) can be run with systems with internal reso- 
nances without additional complexity, something that is often incompatible 
with other formulations [I]. Some comments are given in [T3] regarding the ac- 
cess to higher- dimensional manifolds, but this approach will not be considered 
here. Indeed, resonance frequencies of wind instruments may be proportional 
to each others, in particular if the bore is a cylinder, which is the case for 
the clarinet model, detailed in section [2j The work presented in this paper 
should not be seen as a mere application of known methods to a particular 
example. The NM formulation employed here is very recent, and it is our hope 
that part of the presented work will contribute to the knowledge of how to 
compute the invariant manifolds (see more particularly sections 13.51 and 13 . 6j) . 
In section HJ it is demonstrated how the computation of the NM extending the 
linear modes of the bore, and of their properties, can be used to analyze and 
even predict the model's behavior in terms of transient and steady state, in- 
stantaneous amplitude and frequency, limit cycles and their stability. Results 
are systematically confronted with direct numerical simulations. 



2 THE CLARINET MODEL 



A wide class of musical wind instruments have similar principles of function- 
ing: the player, by blowing inside the instrument destabilizes a valve (a simple 
reed, a double reed or two lips). The acoustic response of the instrument acts 
as a feedback loop which influences the valve behavior. The production of a 
sound corresponds to the self-sustained oscillation of this dynamical system. 
Obviously, in spite of these similarities, the functioning of each class of in- 
struments possesses its own specificities. In this section, basic principles of 
the clarinet functioning are briefly recalled. Simple models are available in the 
literature [T5lfT6lfT71fT8lfT9l . 



2.1 The reed 



The reed is often modelled as a mass/spring/damper oscillator. However, be- 
cause of a resonance frequency (~ 2000Hz) large compared to the first harmon- 
ics of typical playing frequencies, inertia and damping are often neglected [15] . 
This hypothesis leads, considering that reed dynamics is governed by the pres- 
sure difference across the reed, to 

k s (z - Zq) = QOjct - p m outh) (1) 



where z (respectively zq) is the reed position (respectively at rest). The reed 
is closed when z = and opened when z > 0. k s is the reed surface stiffness, 
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Table 1 

Parameter values of the clarinet model Eq. (jlip 



entity 


definition 


value 


unit 


c 


embouchure parameter 


0.35 


1 


7 


blowing pressure parameter 


0.39 


1 


V 


coefficient of viscosity 


0.02 


1 


Yj 


admittance 


1.3r?V2j - 1 


1 


A 


C(3 7 - 1)/(2V7) 




1 


B 


-C(3 7 + l)/(8 7 3/2 ) 




1 


C 


-C(7 + l)/(16 7 5/2 ) 




1 


c 


speed of sound 


340 


m/s 


I 


length of resonator 


0.655 


m 




resonance frequency 


(2j - 1)2ttc/(4Z) 


1/8 




Fig. 1. Scheme of the embouchure of a clarinet. 

Pmouth and pj e t are the pressure deviation in the mouth and under the reed tip, 
respectively. 



2.2 The air flow 



As noted by Hirschberg |20j, in the case of clarinet-like instruments, the control 
of the volume flow by the reed position is due to the existence of a turbulent 
jet. Indeed, a jet is supposed to form in the embouchure (pressure pj et ) after 
the flow separation from the walls, at the end of the (very short) reed channel 
(see Fig. [1]). Neglecting the velocity of air flow in the mouth compared to jet 
velocity fj et , the Bernoulli theorem applied between the mouth and the reed 
channel leads to 

Pmouth = Pjet + ^f™fet ( 2 ) 
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where p is the air density. Since the cross section at the inlet can be expressed 
as the product between the reed opening z and the reed width w r (not visible 
in Fig. [T] since it is transversal to the plane of the figure), Eq. (j2J) can be 
re-written as 



where u is the volume flow across the reed. Combining Eqs. (CEJ) and ([3]) leads 
to the well known expression of the volume flow as a function of the pressure 
difference across the reed 



Since the cross section of the embouchure is large compared to the cross section 
of the reed channel, it can be supposed [2T1 chapter 7] that all the kinetic energy 
of the jet is dissipated through turbulence with no pressure recovery (like in 
the case of a free jet). Therefore, the pressure in the jet is (assuming pressure 
continuity) the acoustic pressure p r imposed by the resonator response to the 
incoming volume flow u. 

Introducing the non-dimensional pressure [22j chapter 6] p = p r / (k s Zo) and 
volume flow u = Z c u/(k s zo), Eq. d3J reads 



u = C(l +P~l)Vl^P (5) 



with ( = ZcUJr^l^f 1 and 7 = p m outh/ '(k s zo). When the reed is closed, i.e. 
1 + p — 7 < 0, the volume flow is zero (u — 0). 

It has been verified [23J that a cubic expansion of Eq. ([5]) leads to a reasonably 
good approximation to the resulting periodic solutions, at least far from the 
complete closing of the reed. Therefore, we assume that the volume flow is 
finally given by 




(3) 



U = W r (z - — (p m outh -Pjet))\/-(Pmouth ~Pjet)- 




(4) 



u = u + Ap + Bp 2 + Cp 3 



(6) 



with u = C(l - t)>/7' A = C 



37-1 
2^ 



B 



c 



37+1 

8 7 3 / 2 



and C 



c 



16 7 5/2 • 



7+1 
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2.3 Acoustics in the instrument 



We consider a cylindrical bore (length /) for the clarinet, and follow Debut [2U 
page 60] . Although the final model will be considered in the time domain, it is 
first written in the frequency domain, where it is simplified before going back 
to the time domain. 

The model retained is the wave equation inside the tube, with a source at 
x = x s to take into account the air flow blown into the instrument, and 
Neumann and Dirichlet boundary conditions at the input and the output of 
the tube, respectively: 



d 2 



(]- + « 

C 



P{w,x) = -]uj^U(uj)5(x s ), Vx G [0,1], 



d x P(x,u) 

P(x,Lu) = 







for x 
for x 




I. 



(7) 



a is a real number representing visco-thermal losses (dispersion is neglected) 
and varies as the square root of the frequency. The case of a clarinet model 
corresponds to a source located at the input of the tube (i.e. x s — ► 0). 

The dimensionless pressure field P(x, ui) is decomposed into the family {/ n }„ e N 
of orthogonal eigenmodes of the air column inside the bore, 

oo 

P(x,C0) = J2fn(x)P n (c0). (8) 
n=l 



In the case of a closed/open cylindrical bore of length / and dispersion ne- 
glected, f n (x) = cosk n (x), where k n = and n is an odd positive integer. 

Modal coordinates P n (cj) are calculated through the projection of Eq. (J7J) 
(with P replaced by Eq. (jSJ), truncated to modes) on each mode f n , leading 
to 

- uo 2 P n {uo) + 2ac ]uP n {u) + {u 2 n - a 2 c 2 ) P n {uo) = ^uU{u) (9) 



where u n = k n c. Several approximations are now made concerning the losses 
coefficient a. Since the damping of each mode is smaH[2H page 70] (ac << u n ), 
the third term of the left hand side of Eq. can be reduced. Moreover, though 
a is a function of the frequency, we consider that a constant value a n can be 
associated to each mode ([2H page 72]). Noting a n ~ ^f, where Y n is the value 
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of the admittance at frequency u n /27[, Eq. (J9]) can now be written in the time 
domain second order ODE: 

jp n {t) + 2Y n jp n (t) + u%p n (t) = jti(t), (10) 



where p n {t) is the inverse Fourier transform of P n {uj). 



2.4 Complete model 



Considering Eqs. (TTU]) and flS]) leads to the dimensionless model, made of N 
second-order ODE, 

p n {t) + 2Y n ±p n (t) + u 2 n p n (t) = 

(11) 



2c N ( N \ \ N 

T U + 25Sft(*) + 3Cf (Eft(*)J JEpM- 



The total pressure at the input end of the instrument is given by the sum of 
all p n (t) (Eq. © with x = 0). 

The intention now is to apply the concept of NM to this system of Eqs. ( TlTT) . 
The aim is to get a reduced representation of the model, which could help in 
analyzing the model behavior, and which could be useful for sound synthesis 
purpose. 



3 NONLINEAR MODES 



It is briefly recalled here how to characterize the NM in the framework of the 
invariant manifold theory using an amplitude-phase transformation according 
to Bellizzi and Bouc [25H3] . 

We consider a system of the form 

MP(t) + F(P(t),P(t)) = (12) 



where P is an n- vector function, M is a non-singular symmetric square N x iV- 
matrix and F is a (sufficiently regular) vector function with dimension N such 
that F(0, 0) = 0. 
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3. 1 Linear modes as a starting point 



In this section, definition and properties of the normal modes are recalled 
when F in Eq. f)12p is a linear function. Variables introduced in this section 
will be then extended to the nonlinear case in further sections. 



3.1.1 Undamped case 

Let us consider the case where F(P, P) = KP, where K is a symmetric square 
iV x iV— matrix. The normal modes are then the N pairs (Q p , *& p ) solutions 
of the eigenvalue problem 

K* p = Mtf p fi* (13) 



and the orthogonality condition and the mass-normalization are written re- 
spectively as 



**M* q =0 Vp^q, (14) 

**M* P =1 Vl<p<n. (15) 

A family of periodic solutions is associated to each normal mode as 

P(t) = vX(<P(t)) (16) 



where 



X(0) = * p cos(0) 

0(t) =n p t + v ■ ( 17 ) 

v = a (constant) 



The amplitude and frequency of the periodic motion are v and Q p /2n respec- 
tively and X is a periodic function with respect to the phase variable 0. 
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3.1.2 Damped case 

If F(P, P) = KP + CP, where K and C are real, square matrices, the eigen- 
value problem to solve is now [27J 



(c 


















°] 






V° 


-M J 



P 



with X p = i] p ± jfip (assuming ^ 0) and *p = (* P T , A* P T ) T With * p 
+ i^fp, the orthonormalization condition can be selected as 

*£ T M*£ + ** T M^ = 1, 

vfM®; = o. 



The family of solutions associated to each mode is now written as 

P(t) = v(t)X(<f>(t)) (18) 

where 

X(0) = & c p cos - *p sin 

<P{t) =Q p t + V ■ (19) 

v(t) = ae r ' pt 



Depending on the sign of r] p , the motion can be damped (r] p < 0), amplified 
(r] p > 0), or periodic (t) p = 0). The amplitude and frequency of the motion are 
v (t) and Q p /2tt respectively. Here also X is a periodic function with respect 
to the phase variable 0. 



3.2 Definition of nonlinear modes 



In the case of a nonlinear function F in Eq. (fT2|) . the linear modal formalism 
recalled in section 13.11 is extended hereafter with the following major differ- 
ences: 

• X is not only a periodic function with respect to the phase variable but is 
also a function of the amplitude v (see Eq. (T2U1) ). 
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• Even in the case of a periodic motion, the amplitude v is a function of time 
defined through a differential equation (see Eq. (12Tb ). The linear damped 
case studied in the above section would correspond to = 77 in Eq. (|2"Th). 

• The pulsation of the motion Q is no more constant, but is a function of the 
amplitude v and the phase (see Eq. (l2~Tb). 

Indeed, we focus on motions (solutions of Eq. (I12p ) where the pressure com- 
ponents and its derivatives (P and P) are related to a single pair of amplitude 
and phase variables {y and (f>) according to 

{l>(t) = v(t)X(v(t),<f>(t)) 
\P(t) = v(t)Y(v(t),cf>(t)) 



where X and Y are N- vector functions, which are 2n -periodic with respect to 
the phase variable <fi and the amplitude and phase variables are governed by 
the two first-order differential equations 

v(t) =v(t)£(v(t),4>(t)) , , 

with ^ . (21) 




In Eq. ([2~Tj) . (the frequency function or frequency modulation function) and 
C, (the damping function or amplitude modulation function) are two scalar 
functions. As established in [13], these two functions can be only chosen as 
even and it-periodic with respect to the phase variable. Furhermore, ip G [0, 27r] 
and a > are two given constants which set the initial conditions of the 
motion. 

If such a family of motions (!20l - (pTI) . parameterized by the variables (a, ip) 
exist, it defines a NM for Eq. (fl2|) . which is characterized by the four functions 
X, Y, Q and f . 



3.3 Some properties of nonlinear modes 



For a given NM, the modal motions are confined to lie on a two-dimensional 
invariant manifold [4,28j in the phase space, defined by the parametric equa- 
tions 

P = aX(a, (p) 

for (a,<p) G R x [0,2tt]. (22) 

P = aY(a,(p) 
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and the modal dynamics on the invariant manifold are given by Eq. (I2ip . 
In terms of signal processing, the function <fi characterizes the instantaneous 
frequency of the modal motion and the function v defines the amplitude mod- 
ulation of the modal motion. 

If the damping function £ = (which means that v(t) = a, Vi), all the modal 
motions (defined by Eqs. f[2"Uj) and fl2~Tj) ) will be periodic. The period is given 
by 

2tt 



r(a) = ' DM)* (23) 







showing that the period is only amplitude dependent. This situation appears 
for autonomous conservative systems |25j. 

Periodic modal motions may also exist if £ ^ or more precisely if 0) 
does not keep a constant sign. Indeed, from Eq. (j2ip . it follows that 

^L = vr{v,(f>) (24) 



where r(v,(j)) = — — - can be viewed as a "generalized damping rate func- 

n(v,(j>) 

tion". Since r is 7r-periodic with respect to the independent variable 0, a 
periodic solution v* (with v*{<p) = v*(<f> + 7r)) may exist for some £ and Q 
(one necessary condition being that r(v,<p) does not keep a constant sign). It 
follows that the associated modal motion 

P(t)=^(0(t))X(t;*(0(t)),0(t)) (25) 



with <f>(t) = £l(v* (<f>(t)) , 4>{t)) and 0(0) = <p, will be T-periodic with period 



2- 



T 



(26) 



From Eq. ( 1241) . the stability analysis of the periodic function v* can be deduced 
using the average principle in the context of perturbation theory [29] from the 
existence of an equilibrium point in the averaged equation 

_ = v < r > (v) (27) 
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where 



More precisely, each equilibrium point v** (defined as < t > (v**) = 0), which 
can be viewed as a constant approximation of a periodic function v*, char- 
acterizes a periodic modal motion (or limit cycle) on the invariant manifold. 
This limit cycle is asymptotically stable if d< ^ v > (v**) < 0. Note that in order to 
analyze the stability in the complete phase space and not only in the invariant 
manifold, the Floquet theory has to be applied [29lfT3] . Finally, the periodic 
modal motion and the associated period can be approximated by, respectively, 

P(t) =v**X(v**,(j)**(t)) (29) 
with 4>**(t) = tt(v**,(j)**(t)) and 0**(O) = <p, and 

2tt 







3-4 Characterization of a nonlinear mode 



Substituting Eq. ( 120]) into Eq. ( !T2l and using Eq. ( 12T]) yields a set of first-order 
nonlinear Partial Differential Equations (PDEs) in the two variables (v,(j>), 

(X + vX v )£ + X (t> n = Y, (31) 
M (Y + vY v ) £ + MY^O + ^-F(vY, vX) = (32) 



where (.)^ and (.)„ denote the partial differentiation with respect to (ft and v, 
respectively. The PDEs fl3T|) - fl32|) are independent of time. 

In order to characterize the four unknown functions (of v and <p) X, Y, Q and 
£, it is necessary to add two scalar constraint equations to fl3TT) and fl32|) (often 
called normalization conditions). Due to the 27r-periodicity with respect to the 
variable 0, the functions X can be expressed as 

X = X oc + X ec + X os + X es (33) 



where (.) oc , (.) ec , (-) os , and (.) es denote the odd cosine, even cosine, odd sine, 
and even sine terms in the corresponding Fourier expansions. We will adopt in 
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this study, without loss of generality, the following scalar constraint equations 
X ocT MX oc + X/ sT MX/ 8 = cos 2 0, (34) 
X ocT MX/ s = 0. (35) 

These constraint equations involve only odd terms of the sine/cosine devel- 
opments due to the assumptions imposed on the two scalar functions Q and 
£ (see Eq. (12~TI) ). Moreover, this choice reduces to usual normalization condi- 
tions for some important special cases [25] including the linear case (see next 
section). 

Finally, a NM of the system (Tl2|) is obtained by solving Eqs. (I3~l~]) - (l35]) for 
the four functions X, Y, Q and £ with initial values given at v = and the 
periodicity properties 

X(v,</>) = X(v,0 + 2vr), 

Sl(v,4>) = Sl(v,-<i>) = n(v,<f> + 'ir), (36) 
£(v,(t>) = £(v,-<f>) = Z(v,<f> + Tr). 

It can be shown [25] that a NM can be defined from each mode of the under- 
lying linear system by selecting it as an initial condition, using the relations 
(1191) . However, contrary to the case for linear systems, an N DOF nonlinear 
mechanical system can possess more that N NM [IT] . 

It is worth noting that depending on the properties of the function F, some 
functions among X os , X es , X oc , X ec can be discarded. For example, 

if F(P, P) = F(P) then X os = 0, X es = , 

if F(P, P) = -F(-P, -P) then X ec = 0, X es = 0. 

3. 5 Numerical solution of the equations describing the manifold 

Eqs. fl3Tl)- fl3"5l) constitute a partial differential algebraic equation (PDAE). 
It is an initial-boundary-value problem, where v acts as a time-like variable, 
and with periodic boundary conditions in the 0-direction. Differential alge- 
braic equations are generally much more difficult to solve than differential 
equations. Firstly, the initial conditions must satisfy not only the algebraic 
constraints, but also a number of compatibility equations depending on the 
index of the equation [30] • For the solutions considered in this work, the NM 
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are the continuations of corresponding linear modes. Consistent initial con- 
ditions of the PDAE are therefore directly obtained from the corresponding 
linearized system, or equivalently by setting v = in Eqs. (I3ip - (l35p . which 
then collapse to an algebraic system of equations. Secondly depending on the 
formulation of the numerical method, there is often a numerical drift in the 
fulfilment of the algebraic constraints as the integration advances. 

The PDAE is solved by the finite difference method. The unknowns are dis- 
cretized in the ^-direction, after which the occurring derivatives with respect 
to are approximated using finite difference approximations. This so-called 
method of lines [31], where the v -derivatives are still left on their continu- 
ous form, leaves us with a usual differential algebraic equation (DAE), which 
can be solved with a suitable numerical integration scheme in the f-direction. 
The backward differentiation formulae (BDF) are a wide class of methods for 
DAEs, of which the implicit Euler scheme is the most well-known represen- 
tative. The approach differs from those used by Pesheck et al. [32] or Bellizzi 
and Bouc [13] , where use is made of Galerkin methods based on trigonometric 
terms, and in the latter case polynomial terms in the f-direction. Although 
elegant, the Galerkin treatment becomes prohibitively complex as the number 
of expansion terms increases. A low order, implicit scheme for step-wise ad- 
vancement in the w-direction is also better adapted for capturing variations, or 
even irregularities, in the f-direction of the solution, which do not readily lend 
themselves to an accurate description with a polynomial basis. A step-wise in- 
tegration is also consistent with the initial value character of the equation 
in the sense that at each v, the solution depends on the earlier, but not on 
subsequent states. 

Let Z(v,<p) be anyone of the unknowns £(i>, 0), fl(v, <f>), X l (v,4>), Y l (v,<f)) for 
i = 1, . . . , n and X(v, (p) = {X x {v, (j)),X 2 {v, 0), . . . , X n (v, <p)) T . The approxi- 
mation U is defined through the discretization 



Z(vj, <f> k ) « U j)k , Vj = jh v , j = 0, 1, . . . , 

(p k = kh (f) , k = 0,1, - 1, = 

In the sequel, the respective discrete approximations Uj t k are denoted X^ k , 
Yf k etc. where the meaning is clear. Since the problem is periodic in the <fi- 
direction, and can be expected to have a smooth solution in this direction, it 
is natural to use a pseudo-spectral [26J approximation of the <9/<9</>-terms. The 
approximation can be interpreted as a usual finite difference approximation 
where the number of points, and hence the order of accuracy, is a function of 
the step-size h$ in the direction of differentiation. The implementation involves 
manipulations in the Fourier space and relies on the fast Fourier transform 
applied to the grid data. Accordingly, a term dZ(v k ,(j)k)/d(j) is approximated 
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by the difference scheme 



AV2 



9Z(vj,(j)k) _ n „ A v , T j 

V m=-JV <A /2 



where .D^ is denoted the pseudo-spectral differential operator, and the coeffi- 
cients d m depend on the discretization. The second index of U reflects the fact 
that data from a grid — periodic in the ^-direction — is used in a wrap-around 
fashion. For the appropriate choice of d m , the approximation has spectral con- 
vergence rate, meaning that the error decreases faster than any polynomial in 

v 

Leaving the system of equations (|3T1) - (!35|) stated on its implicit form, and by 
approximating the appearing derivatives with the pseudo-spectral scheme in 
the 0-direction and a backward difference in the w-direction, the implicit Euler 
approximation of the PDAE at hand is given by 

x j+i,k + v j+i 3+1 i v — — + D ^ x j+i,k^j+i,k = Y j+i,ki 

i = 1, ... k = 0, . . . , N<f, - 1 

2^=i Mu \ Y j+l,k + v j+i JT V Sj+l,k + ^*j+i,fc"j+l,fc \ H 

i = 1, ... ,n, fc = 0, . . . - 1 

X°f 1>fe [M]X- 1>fc + D^f^MD^X™^ = cos 2 fc G {LI} 
X°f ljfe [M]^X^ 1>fc = 0, G {LI}. 



The decomposition in even and odd cosine and sine parts is done with the 
aid of the discrete Fourier transform (DFT), followed by a selection of the 
appropriate Fourier components, and finally an inverse DFT. Due to the sym- 
metric properties of the unknowns, not all the equations k = 0, . . . , — 1 
are linearly independent. This reduces the number of equations, but since it 
is known beforehand that Q and £ can be expanded in only even cosine terms, 
it is possible to reduce the number of unknowns correspondingly. Accordingly, 
the set {LI} is chosen so that only linearly independent equations are retained, 
and new variables £ ec and Q ec , containing only the odd cosine components of 
£ and Q, are introduced as unknowns in the numerical solution of Eq. (1371) . 

The solution on each new f-level j + 1 is obtained by solving Eq. (13T1) for all 
Aj +l fc , Yj +lk , and This non-linear system of equations is solved 

numerically using the Newton method, where the starting solution in each 
step is obtained as a first order extrapolation of the solution at levels j — 1 
and j. 
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3. 6 Computation of time dependent solutions from the manifolds 



Once the surfaces X(t>,0), Y(t>,0), £(t>,0) and fl(v,(p) have been computed, 
the time evolutions v (t) and <fi{t) can be computed numerically solving Eq. (12Tj) . 
Since numerical approximations of £ and Q are known only for certain discrete 
values Vi and 4>u a two-dimensional interpolation procedure is used, employing 
trigonometric interpolation in the ^-direction and quadratic interpolation in 
the f-direction. For solutions representing a limit cycle, it is necessary to solve 
Eqs. ( JHT|) - (j55|) on an interval [0, f max ], where t> max is large enough so that the 
computed region of the invariant manifold contains the limit cycle. In practice, 
this implies a value of t> max slightly above the amplitude of the limit cycle v** 
as defined in section 13.31 An estimate of when t> max is reached can be obtained 
by keeping track of the mean damping function Eq. (1281) . which is a measure 
of the average energy dissipated or supplied to the system over one period. 
It is zero at v — v**. The physical variables in phase space are finally given 
by Eq. fl20|) . where once again the two-dimensional interpolation procedure is 
employed. 



4 Nonlinear modes for the clarinet model 



The clarinet model described in section [2] is considered for a case where N = 3. 
As the excitation of modes four and onwards is fairly weak for the chosen 
blowing pressure, only the most prominent three modes are treated for clarity. 
The method presented in section [3] is applied to find the NM of the system, 
with P = [pi,P2,P3\ T ■ With Eq. (TTTj) linearized and written on first-order 
form, it is possible to determine the value of the blowing parameter 7 where 
an equilibrium point loses its stability and a self-sustained oscillation can 
appear if a Hopf-bifurcation occurs. The motion becomes oscillatory when the 
corresponding eigenvalue of the linearized system matrix crosses the imaginary 
axis. The model parameters (Tab. [1]) are chosen so as to correspond to a mezzo 
forte playing condition. The first vibrational mode becomes linearly unstable 
at 7 = 0.363, see Fig. [2J For 7 = 0.386, also the second mode becomes linearly 
unstable. Thus, the chosen blowing pressure (7 = 0.39) is just strong enough 
to render also the second mode linearly unstable. It means that for the chosen 
value of 7, playing in the first or second register might be possible, depending 
on their stability and on the initial conditions. This will be investigated in the 
following. 
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Fig. 2. Real parts of the eigenvalues of the system matrix of a linearized, first order 
version of Eq. (jllh as a function of 7 (with parameter values given in Tab. [TJ). 

4-1 First mode 



In a first step, the NM described by X and Y, as well as £ and fi, is computed 
solving Eqs. fl3Tl) - fl35l) with the method described section 13.51 and by choosing 
the first linear mode of the linearized system (which differs from the first 
linear mode of the resonator) as the initial condition at v — 0. The shapes of 
the surfaces X, £ and Q are shown in Fig. [3] for a case with 31 discretization 
points in the ^-direction (effectively resolving the first 15 Fourier terms), and 
50 discretization points on the interval v e [0,0.41]. As expected, the shape 
of X\ starts out from a purely harmonic variation with respect to (j) at v = 0. 
The shape then changes only slightly as the amplitude grows. Components X 2 
through X3 are small at v — 0, but then change dramatically as the amplitude 
grows. For low amplitudes, £ is constant and positive, which is characteristic 
for an unstable motion. As the amplitude grows, the shape becomes more 
complicated attaining both positive and negative values. The function < r > 
(Eq. ([HD ) starts out at a positive value, and crosses the v — axis at v** = 0.4046 
with a negative derivative. This is a necessary condition for a self-sustained 
oscillation (see section [3731) . 

In a second step, a time evolution v(t) and <p(t) is calculated from an arbitrary, 
small-amplitude initial condition (v (0) = 0.1, 0(0) = 0) numerically solving 
Eq. ( I2TI) . as described in section [3751 According to the obtained results, the 
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time evolutions of the components of P(t) and P(t) can finally be computed 
with Eq. ( 1201) . These results are shown in Fig. HHH1 Instead of <p(t), the in- 
stantaneous frequency <fi(t) is shown. As expected, for small values of v, <p 
oscillates around 815.55 rad s^ 1 which is the value of the imaginary part of 
the first eigenvalue of the linearized system. It is noticeable that there is a 
slight difference between this frequency and the corresponding resonance fre- 
quency of the 1/4— wavelength resonator (27rc/4/ = 815.37 rad s -1 ). For larger 
values of v, the amplitude of the oscillation in <fi is increasing. As can be seen 
in Fig. [J], the amplitude of the first component p\ grows quickly initially, but 
then stabilizes as the limit cycle is approached. The same phenomenon is visi- 
ble for the higher components, but they show a much stronger relative growth. 
This is a typical feature for wind instruments, where the small amplitude os- 
cillations are nearly sinusoidal. As the amplitude grows, nonlinear effects add 
increasingly to the timbre by successive enrichment of the harmonic content 
of the signal. The envelopes of the components p n are defined by the form of 
v together with the evolution of the components of X versus v. 

In order to check the validity of the solutions, a reference solution (with initial 
conditions given by Eq. (1201) at t = 0) was computed by direct solving of 
Eq. (fTTj) using a Runge-Kutta-Fehlberg solver (ode45 in Matlab) with a small 
tolerance. The difference between the solution obtained from the MN approach 
and the reference solution is presented as the error in Fig. [5j The error is very 
small during the transient phase, but then starts to grow slowly in time mainly 
due to the error in the frequency. A longer simulation demonstrates (not shown 
in figure) that the error envelope grows for some time approximately linearly 
with time, consistent with a slowly growing phase lag. (The same kind of error 
growth would be prevalent for any method with numerical dispersion, however 
small the error in the frequency, albeit at a different rate.) 

It is interesting to compare the limit cycle computed from the NM, with that 
obtained from the reference solution. The comparison eliminates any accumu- 
lated phase errors, and gives another more direct estimate of the error. Fig. [6] 
shows the limit cycles superimposed. The difference is hardly distinguishable 
in the plot. An approximation of the limit cycle computed from Eq. ff29l) . with 
the constant value v** = 0.4046 (zero of < r >, see bottom left of Fig. [3D, is 
also shown in the figure. Evidently, a good approximation of the limit cycle 
can also be obtained. In order to get an estimate A p of the amplitude of the 
pressure signal p = £f =1 Pi in the steady state, surfaces Xi must be taken into 
account through Eq. (129!) 




The limit cycle frequency given by Eq. ([301) is 815.2 rad s l , which can be 
compared to the linear resonance frequency of the 1/4 wave resonator that 
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(a) (b) 




v 



Fig. 3. First nonlinear mode of the clarinet: (a) Damping function £, (b), (d) and 
(f) Surfaces for the three components of X, (c) Instantaneous frequency f2, (e) The 
scalar function < r >, the zero of which indicates an estimate of the amplitude of 
the limit cycle. 

is 815.37 rad s^ 1 , and to the eigenfrequency of the linearized system that is 
815.55 rad s _1 . Thus a modal motion according to the first NM corresponds 
to a clarinet sounding in its first register. We have seen that the instantaneous 
amplitude and frequency are oscillating, even in the steady state. This is not 
in contradiction to the fact that the corresponding regime in terms of time 
evolution of pi is periodic. This can be checked from Fig. [6]where the dynamics 
of the steady state corresponds to a limit cycle in the configuration space. 

4-2 Second mode 

A similar investigation for the second nonlinear mode of the clarinet is pre- 
sented in Fig. [TrfTUl The initial conditions for the computation of the nonlinear 
mode have been changed to the second linear mode of the linearized system. 

We see in Fig.[7]that the second component X2 is now dominating over X\ and 
X3. The value of £ shows a variation that initially resembles that of the first 
mode, with positive as well as negative values. The amplitude v** = 0.1035 of 
the limit cycle is again found as the zero of < r >, which is smaller than for 
the first mode. Also for the second mode, the derivative of < r > is negative at 
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(a) v(t) 




0.1 0.2 0.3 0.4 0.5 



t(s) 

Fig. 4. First nonlinear mode: (a) Time evolution of v, (b) Time evolution of (f>. 
Barely discernible is a fine ripple in v with the same period as the limit cycle. 
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Fig. 5. First nonlinear mode: (a), (c) and (e) Time evolution of pi, p2 and 
characterizing a modal motion according to the NM approach, (b), (d) and (f) 
Errors computed with a direct simulation of Eq. (jlip and the same initial conditions 
as reference. 
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0.4 -0.1 



Fig. 6. First nonlinear mode: The limit cycle in the configuration space (pi,P2,P3)- 

- Limit cycle computed from the nonlinear mode, limit cycle according to 

direct simulation, limit cycle computed from v** . 

v = v**, showing that the limit cycle is stable on the invariant manifold. This 
means that the limit cycle is stable with respect to a subspace of disturbances, 
but not necessarily to any disturbance. 

Fig. [8] shows the time evolution of v and <fi calculated from an arbitrary, 
small-amplitude initial condition (v (0) = 0.01, 0(0) = 0) solving numerically 
Eq. ( I2ip . It is noticeable that ripples on v and <fi are much weaker than in 
Fig. HI This is related to the smoothness of the surfaces £ and Q (see Fig. [7|), 
which remain much more regular than for the first mode, even past the limit 
cycle amplitude. 

Fig. [9] shows the limit cycle of NM 2 on the invariant manifold represented 
by the three components pi~ps according to Eq. (1291) . The smoothness of the 
invariant manifold is linked with the smoothness of the surfaces of X, which 
are again much smoother even for high amplitudes, than is the case for the 
first mode. The invariant manifold for the first mode would be much more 
difficult to visualize due to its intricate folds and intersections with itself. The 
shown surface is really a projection of a six- dimensional manifold (including 
also P1-P3) that does not intersect with itself. 

A direct numerical simulation (solving Eq. pip ) indicates that the limit cycle 
for the second nonlinear mode may be unstable in the phase space for the 
chosen parameter values. Indeed, a direct numerical simulation started from 
initial conditions on the manifold representing the second nonlinear mode, will 
eventually jump out of the manifold due to round-off and truncation errors. 
This is seen from the curves in the center column of Fig. [TO], The time evolution 
is initially the same as for the solution computed from the NM approach 
(left column), but after about 1.5 s, the first component starts to become 
notably excited and the growth of the second component is discontinued. The 
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Fig. 7. Second nonlinear mode of the clarinet: (a) Damping function £, (b), (d) and 
(f) Surfaces for the three components of X, (c) Instantaneous frequency 0, (e) The 
scalar function < r >. 

limit cycle for the second nonlinear mode is never reached, and the oscillation 
eventually converges to the limit cycle of the first mode. 

To analyze the stability of the limit cycle, Floquet Theory has to be applied 
|29j . Starting from the approximated limit cycle given by Eq. (I29I) . the mon- 
odromy matrix is computed solving the 27r-periodic variational linear differ- 
ential system associated to Eq. (TTT1) over one period, using the six canonical 
basis vectors as initial conditions successively (see details in [13]). The com- 
putations show that the periodic orbit (approximated by Eq. fl29|) ) associated 
with the second NM is unstable in the phase space (two complex conjugate 
multipliers are outside the unit circle). 

Thus, for this clarinet model, the second register appears unstable with the 
chosen parameter values, which is often experienced by beginners on real in- 
struments. 



4-3 Third mode 



Calculations for the third NM are presented in Fig. [TTHT21 The initial condi- 
tions for the computation of the NM have been changed to the third linear 
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Fig. 8. Second nonlinear mode: (a) Time evolution of v, (b) Time evolution of (f>. 




x 10 3 



Fig. 9. Second nonlinear mode: The limit cycle in the configuration space (pi,P2,P3) 
lying on the invariant manifold. Computations are carried out from the NM ap- 
proach. 

mode of the linearized system. 

We see, Fig. [TH that the third component X3 is now dominating over X\ and 
X 2 . The surface £, unlike for the first and the second modes, starts out at 
a negative value for v — 0. As v increases, the surface becomes increasingly 
oscillatory, but at no point is the mean value positive. In terms of the mean 
damping function, < r > is strictly negative. As a consequence, for the cho- 
sen blowing parameter 7, a solution started at any point on the third mode 
invariant manifold converges to the equilibrium point which is stable. There 
is no limit cycle in the invariant manifold associated to the third NM. This 
is exemplified, in the time domain on v and <p with small-amplitude initial 
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(d) p^(NM) (c) p ^ (direct simulation) (f) \difference for p^\ 
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t (s) t (s) t (s) 

Fig. 10. Second nonlinear mode: (a), (d) and (g): Time evolution of pi, p2 and p% 
characterizing a modal motion computed from the NM approach, (b), (e) and (h) 
Time evolution computed by direct simulation of Eq. (11 If) using the same initial con- 
ditions, (c), (f) and (i) Differences between solutions computed in the two different 
ways. 

condition v(0) = 0.1 and 0(0) = (see Fig. [12]). 



4-4 Discussion 

Some problems with divergence of the solution in the computation of the 
NM have been observed for certain values of v. Numerical experiments indi- 
cate that the value t>di v , where the instability occurs, converges to a certain 
value as the step size h v in the w-direction decreases. In the case of a model 
with one single degree of freedom, it is possible to formulate a version of 
Eqs. (I3l]) -(j32~ |) where the algebraic constraints Eqs. ( |34|) -(!35l) are eliminated 
by assuming ~K(v(t),4>(t)) = cos(0(i)) in Eq. ( 1201) . The divergence persists 
also in that case, which indicates that the problem is not a consequence of the 
differential-algebraic structure of the original problem. The findings suggest 
some sort of ill-posedness in the continuous equation, rather than a numerical 
problem, but the exact nature of this is yet to be investigated in detail. For 
higher blowing levels, the point of instability is reached before the amplitude 
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(a) ' ...... (b) 




v 



Fig. 11. Third nonlinear mode of the clarinet: (a) Damping function £, (b), (d) and 
(f) Surfaces for the three components of X, (c) Instantaneous frequency Q, (e) The 
scalar function < r >. 




Fig. 12. Third nonlinear mode: (a) Time evolution of v, (b) Time evolution of 4>. 

of stable oscillation is reached, thus limiting the applicability to medium am- 
plitudes. However, the validity of the clarinet model itself is limited to medium 
amplitudes (i.e. non beating reed). 

The clarinet model with up to 8 degrees of freedom has been investigated, with 
results that are in accordance with the case N = 3. There is in principle no 
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limit of N, but the computational time increases with the size of the system. 



5 Conclusion 

We have seen that the suggested method for computing NM motions of dy- 
namical systems, using amplitude and phase as master variables, is capable of 
dealing with self oscillating systems with internal resonances. The numerical 
method presented in this paper offers a flexible way to handle the problem 
of numerical refinement for increased accuracy, and allows for the treatment 
of larger systems with many degrees of freedom - unlimited in principle, but 
limited by memory demands and execution time considerations in practice. A 
bottleneck in the computations is the solution of Eqs. (|3TI) - (l35l) . The implicit 
method requires the solution of a nonlinear system of equations, whose size 
grows with N and N^, but the use of an analytical Jacobian matrix speeds up 
the computation. Although the calculation of the manifolds is unwieldy, it is 
a pre-processing step after which the computational work is independent of, 
or grows only linearly (the work of forming the components pj) with N. The 
surfaces X, Y, £ and O are intricate for amplitudes in excess of the limit cycle 
amplitude v**, but for smaller v, their smoothness allows for a more compact 
reduced order representation, e.g. each surface might be represented by one 
single, or a piece-wise patchwork, of approximating functions with a limited 
number of parameters. This could greatly reduce storage demands and the 
need for interpolation. 

This approach has many valuable by-products. For example, after having 
solved Eq. (l2"Tj) . which is an ODE with two unknowns, one has immediate 
access to the instantaneous amplitude and frequency without the need to 
solve the whole system (TTTT) . Moreover, one has also direct access to the solu- 
tion components Pi,...,Pn an d Pi, ■ ■ ■ ,Pn for any t. We have also seen that 
is possible to compute unstable dynamics, including transients and the steady 
state. The limit cycles, finally, can be estimated without any computation 
whatsoever in the time domain. However, the manifold (and the dynamics) 
are calculated for constant parameters 7 and (. Therefore, a consequence high- 
lighted by this study is that a direct application of the approach for sound 
synthesis is therefore not obvious. On the contrary, the approach is predomi- 
nately adapted for analyzing model dynamics. 



Acknowledgments 

The authors want to thank Jean Kergomard for fruitful discussions. 



26 



The study presented in this paper was lead with the support of the French Na- 
tional Research Agency ANR winthin the CONSONNES project. The first au- 
thor was financed by the French Ministry of Research through a post-doctoral 
grant. 



References 

[I] H. L. F. Helmholtz, "On the sensation of tone" (Masson) (1877), reprinted by 
Dover, from the first edition in 1954. 

[2] N. H. Fletcher and T. D. Rossing, The Physics of musical instruments (Springer 
Verlag) (1991). 

[3] R. Rosenberg, "The normal modes of nonlinear n-degree-of-freedom systems" , 
Journal of Applied Mechanics 29, 7-14 (1962). 

[4] S. Shaw and C. Pierre, "Non-linear normal modes and invariant manifolds", 
Journal of Sound and Vibration 150(1), 170-173 (1991). 

[5] S. Shaw and C. Pierre, "Normal modes for nonlinear vibratory systems", 
Journal of Sound and Vibration 164(1), 85-124 (1993). 

[6] F. Pellicano and F. Mastroddi, "Applicability conditions of a nonlinear 
superposition technique", Journal of Sound and Vibration 200(1), 3-14 (1997). 

[7] D. Jiang, C. Pierre, and S. Shaw, "Large-amplitude non-linear normal modes of 
piecwise linear systems", Journal of Sound and Vibration 272(3-5), 869-891 
(2004). 

[8] L. Jezequel and C. Lamarque, "Analysis of nonlinear dynamical systems by 
normal form theory", Journal of Sound and Vibration 149(3), 429-459 (1991). 

[9] A. Nayfeh, Method of normal forms (John Wiley Sons, New York) (1993). 

[10] C. Touze, O. Thomas, and A. Chaigne, "Hardening/softening behaviour in non- 
linear oscillations of structural systems using non-linear normal modes" , Journal 
of Sound and Vibration 273(1-2), 77-101 (2004). 

[II] A. Vakakis, "Nonlinear normal modes (nnms) and their applications in vibration 
theory: an overview", Mechanical Systems and Signal Processing 11(1), 3-22 
(1997). 

[12] A. Vakakis, L. Manevitch, Y. Mikhlin, V. Pilipchuk, and A. Zevin, Normal 
modes and localization in nonlinear systems (Wiley Interscience, New York) 
(1996). 

[13] S. Bellizzi and R. Bouc, "An amplitude phase formulation for nonlinear modes 
and limit cycles through invariant manifolds" , Journal of Sound and Vibration 
300(3-5), 896-915 (2007). 



27 



[14] S. Bellizzi and R. Bouc, "Analysis of multi-degree of freedom strongly non-linear 
mechanical systems with random input. Part I: Non-linear modes and stochastic 
averaging", Probabilistic Engineering Mechanics 14, 229-244 (1999). 

[15] J. Backus, "Small vibration theory of the clarinet", J. Acoust. Soc. Amer. 35, 
305 (1963). 

[16] C. J. Nederveen, Acoustical aspects of woodwind instruments (Fritz Knuf pub., 
Amsterdam) (1969). 

[17] W. E. Worman, "Self-sustained nonlinear oscillations of medium amplitude in 
clarinet-like systems", Ph.D. thesis, Case Western Reserve University (1971), 
Ann Arbor University Microfilms (ref. 71-22869). 

[18] T. A. Wilson and G. S. Beavers, "Operating modes of the clarinet", J. Acoust. 
Soc. Amer. 56, 653 (1974). 

[19] R. T. Schumacher, "Ab initio calculations of the oscillations of a clarinet", 
Acustica 48, 71-85 (1981). 

[20] A. Hirschberg, J. Gilbert, A. P. J. Wijnands, and A. M. C. Valkering, "Musical 
aero-acoustics of the clarinet" , Journal de Physique IV 4:C5-559:C5-568 (1994), 
colloque C5 supplement au Journal de Physique III. 

[21] A. Hirschberg, J. Kergomard, and G. Weinreich, eds., Mechanics of Musical 
Instruments, chapter Aero-acoustics of wind instruments, by A. Hirschberg 
(Springer Verlag) (1995). 

[22] A. Hirschberg, J. Kergomard, and G. Weinreich, eds., Mechanics of 
Musical Instruments, chapter Elementary considerations on reed-instrument 
oscillations, by J. Kergomard (Springer Verlag) (1995). 

[23] C. Fritz, S. Farner, and J. Kergomard, "Some aspects of the harmonic balance 
method applied to the clarinet", Applied acoustics 65, 1155-1180 (2004). 

[24] V. Debut, "Deux etudes d'un instrument de musique de type clarinette : 
analyse des frequences propres du resonateur et calcul des auto-oscillations 
par decomposition modale (Two studies of a clarinet-like musical instrument : 
analysis of the eigen frequencies and calculation of the self-sustained oscillations 
by modal decomposition)" , Ph.D. thesis, Aix-Marseille 2 university (2004), URL 
|http : //tel . ccsd. cnrs . f r/tel-00008711| . 

[25] S. Bellizzi and R. Bouc, "A new formulation for the existence and calculation 
of non-linear normal modes", Journal of Sound and Vibration 287(3), 545-569 
(2005). 

[26] D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods: Theory 
and Applications (Society for Industrial and Applied Mathematics) (1977). 

[27] L. Meirovitch, Computational Methods in Structural Dynamics (Sijthoff 
Noordhoff) (1980). 



28 



[28] R. Arquier, S. Bellizzi, R. Bouc, and B. Cochelin, "Two methods for the 
computation of nonlinear modes of vibrating systems at large amplitudes" 
Computer and Structures 84, 1565-1576 (2006). 

[29] J. Hale, Ordinary Differential Equations (Wiley- Interscience, New York) (1969). 

[30] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff and 
Differential- Algebraic Problems (Springer Verlag) (1996). 

[31] B. Gustafsson, H. O. Kreiss, and J. Oliger, Time dependent problems and 
difference methods (John Wiley & Sons, New York) (1995). 

[32] E. Pesheck, S. Shaw, and C. Pierre, "A new Galerkin-based approach for 
accurate nonlinear normal modes through invariant manifolds", Journal of 
Sound and Vibration 249(5), 971-993 (2002). 



29 



